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Abstract 



f^ ' We present a spectral method for parabolic partial differential equa- 

t~^ , tions with zero Dirichlet boundary conditions. The region SI for the prob- 

^O ' lem is assumed to be simply-connected and bounded, and its boundary is 

»«fs ' assumed to be a smooth surface. An error analysis is given, showing that 

(— ^ , spectral convergence is obtained for sufficiently smooth solution functions. 

fSJ ' Numerical examples are given in both R^ and R^. 

1 INTRODUCTION 

X 

$H ■ Consider solving the parabolic partial differential equation 



^- E ^(«m(M,.(M))^)+/(M,.(M)), (1) 

for s e fi C R'', < t < r. The solution u is subject to the Dirichlet boundary 
condition 

u{s,t) = 0, sedn, 0<t<T (2) 

and to the initial condition 

u(s,0) = Uo(s), sen. (3) 



The region 17 is open, bounded, and simply connected in R'^ for some d > 2, 
and the boundary dQ is assumed to be several times continuously difFeren- 
tiable. This paper presents a spectral method for solving this problem. The 
functions aij {s,t,z) and f{s,t,z) are assumed to be continuous for {s,t,z) G 
H X [0,T] X R. Additional assumptions are given later in the paper. These as- 
sumptions are stronger than needed for the results we obtain, but they simplify 
the presentation. In addition, we assume that there is a unique solution u to 
the problem (II])-©. For an introduction to the theory of nonlinear parabolic 
problems using variational methods, see [26l Chap. 30]. 

We transform the above problem to one over the unit ball B^ in M'^, and 
then we use Galerkin's method with a suitably chosen polynomial basis to ap- 
proximate the solution u. This is similar in spirit to earlier work in [2], [5], [7]. 
This approach reduces the problem to the solution of an inital value problem 
for a system of ordinary differential equations, for which there is much excel- 
lent software. The convergence analysis of the paper depends on the landmark 
paper of Douglas and Dupont [T3] . The methods of this paper also extend to 
having the functions atj and / depend on the first derivatives du/dsj, although 
this is not considered here. For related books on spectral methods for partial 
differential equations, see [10]- [12], [16], [H], [22], [^■ 

The spectral method is presented and analyzed in §2, implementation issues 
are discussed in §3, and numerical examples in K.^ and M.^ are given in §4. 

2 A spectral method 

Wc transform the problem (II|)-([3]) to one over the unit ball B^, and then we ap- 
ply Galerkin's method using multivariate polynomials as approximations of the 
solution. To transform a problem defined on J7 to an equivalent problem defined 
on Bcj, we review some ideas from [5] and [7], modifying them as appropriate 
for this paper. 

Assume the existence of a function 

$ : Id i4 n (4) 

onto 

with $ a twice-differentiable mapping, and let ^ = $~^ : il — > M^. For 

onto 

v{x) = w ($ (x)) , a; e Bd C K'' (5) 

and conversely, 

vis) = v{^'{s)), seUcW^. (6) 

Assuming v E H^ (Jl), we can show 

V^v{x) = J{xfVsv{s), .s = $(x) 

with J [x) the Jacobian matrix for $ over the unit ball B^, 

d 

run: t If 1 

J{x) = (D^) {x) 



dxj 



X e Bd. (7) 



To use our method for problems over a region J7, it is necessary to know explicitly 
the functions $ and J. We assume 

det J{x) 7^ 0, x€ Irf. (8) 

Similarly, 

\/sv{s) = K{s)^V^v{x), X = *(s) 

with K{s) the Jacobian matrix for ^ over 51. By differentiating the identity 

*($ (x)) = X, X eMd 

we obtain 

if ($(a;)) = J(x)"^ 

Assumptions about the differentiability of v (x) can be related back to assump- 
tions on the differentiability of ti(s) and ^{x). 

Lemma 1 // $ e C" (Sd) and v £ C^ (fl) , then v e C^ (Bd) with q = 
minjfc, m}. 

Proof. A proof is straightforward using ([5]) . ■ 

A converse statement can be made as regards v, v, and \E' in ([5]). 

Often a mapping (p is given from S''"^ onto dfl, and it will not be clear as 
to how to extend the mapping to $ satisfying Q and ([8]). This is explored in 
[6] with several methods given for constructing $. 

To obtain a space for approximating the solution u of our problem, we pro- 
ceed as follows. Denote by n„ the space of polynomials in d variables that are 
of degree < n: p G n„ if it has the form 

\i\<n 

with i a multi-integer, i = (zi, . . . , id), and |i| — ii + - ■ ■ + id- Our approximation 
space with respect to B^ is 

Xn = { (i - l^f) Pix) I P e n„} C iji {Md) (9) 

With respect to Q, the approximating subspace is 

Xn = {i' is) - ^(* is)) : V^ e a;} c ijgi (n) (10) 

Let Nn = dim <%"„ = dimi;^ = dimH,,. For d ^ 2, N„ = (n + I) {n + 2) /2. 



2.1 The approximation 

We reformulate the parabolic problem ((H)-® as a variational problem. Multi- 
ply ([T|) by an arbitrarily chosen v E Hq (fl) and perform integration by parts, 
obtaining 

du(-,t) \ v-^ f , , .. du(s,t) dv(s,t) , 

+ {fi;t,ui;t)),v), VeH^in), i>0. 

In this equation, (•, •) denotes the usual inner product for L^ (fl) . Equation 
(jlip . together with (|3]), is used to develop our approximation method. 
We look for a solution of the form 

Un{s,t)^Y,akit)^Pk{s) (12) 

fc=l 

with {ipi, . . . , iPn} a basis of <%"„. The coefficients {ai, . . . , aN„} generally will 
vary with n, but we omit the explicit dependence to simplify notation. Substi- 
tute this Un into (fTTj) and let v run through the basis elements ipi. This results 
in the following system: 



k=l 

= -X^a/c (t) Y^ / aij s, t, ^Qffe (t) V'fc (s) I 
fc=i ij=i"'f^ \ fc=i / 



diJk{s,t) dipt{s,t) 
dsi dsj 



I • ■ • I -"m 



(13) 

This is a system of ordinary differential equations for the coefficients a^, for 
k = 1, . . . , Nn- For the initial conditions, calculate 

uo{s) Ri uo,n{s) = Y^a],'tpk{s) (14) 

fe=l 

by some means, and then use 

afe(0) = 4"\ fc = l,...,7V„. (15) 

The implementation of p^ - p^ is discussed in ^ 

2.2 Convergence analysis 

Our error analysis of (|T2|) - p5|) is based on Douglas and Dupont [13l Thm. 7.1]; 
and as in that paper, we assume the functions {atj} and / satisfy a number of 
properties. 



Al As stated earlier, we assume the functions Uij (s,t, z) and f (s,t,z) are 
continuous for {s,t,z) S fi x [0,T] x R. Moreover, assume 

\f{s,t,r)-f{s,t,p)\<K\r~p\, 

for all (s, t, r) , (s, t,p) eUx [0, T] x R, and 

|a»j- {s, t, r) - a,,j (s, t,p)\ < K \r ~ p\ 

for all (s,i,r),(s,i,p) e H x [0,r] x R, 1 < i, j < d. 

A2 We assume that the matrix A{s,t,z) = [ciij {s,t, z)]^ -^ is symmetric, 
positive definite, and has a spectrum that is bounded above and below by 
positive constants 771 and 772, uniformly so for {s,t,z) Gftx [0,T] x R. 

Theorem 2 (Douglas and Dupont) Assume the functions Oij (s, t, z) and f (s, t, z) 
satisfy the conditions A1-A2. Let u be the solution of (QP-jO^ '^'^'^ assume it is 
continuously differentiable over Vt x [0,T]. Let Un be the solution of I J^ |] - JJ5)) . 
Then there are positive constants 7 and C for which 

llu-UnlliaxL- + ^ll" - "nil ffi xL^ 

< C {\\uq ~ U0,n|li2 + \\u - w||i2xL~ (16) 

for any w of the form given on the right side of (J 
The norms used in (|16p are given by 



IMIl^xl- = sup \\v{-,t) \\L2{n) 

0<t<T 

||w||l2xL2 = ||w||L2(nx[0,T]) 

IMl^xL^ = / \\'"i-^t)\\H},{n)dt 
^ 

The assumptions of the theorem imply the assumptions used in |131 Thm. 7.1], 
and the conclusion follows from the cited paper. 

To apply this theorem, we need bounds on the norms given in (I16|) for u — w. 
To obtain these, we use the following approximation theoretic result that follows 
from Ragozin pO] . 

Lemma 3 Assume that g (x, t) , dg {x, t) /dt are k times continously differen- 
tiable with respect to x £ B^;, for some k > and < t < T. Further, 
assume that all such k^^-order derivatives satisfy a Holder condition with ex- 
ponent 7 e (0, 1] and with respect to x £ B^, 

\h{x,t) -h{y,t)\ < Ck,',{g)\x-y\'^ , 
dh{x,t) dh{y,t) 



dt dt 



< ckn (ff) \x - y\ 



uniformly for x,y G M^ and < t < T, where h denotes a generic k*^-order 
derivative of g with respect to x ^Md- The quantity Ck,-y (g) is called the Holder 
constant. Let {cpi, . . . , lpn} denote a basis of Tin- Then for each degree n > 1, 
there exists 

9n (x, t) =^/3k {t) fk (x) 
fe=l 
which satisfies 



max max 



dg (x, t) dgn [x, t) 



dt dt 

for some constant bk.-y > that is independent of g. 






Proof. This result can be obtained by a careful examination of the proof of 
Ragozin [201 Thm. 3.4]. A similar argument for approximation of a parameter- 
ized family g {x,t) over the unit sphere S'^"^ is given in [S]. The present result 
over Brf follows by combining that of [HI §4.2.5] over S'^ with the argument of 
Ragozin over B^. ■ 

Next, we must look at the approximation of the solution u{x,t) by means 
of polynomials of the form given on the right side of (IT^ . To do this, we use a 
trick from [2, (9)-(15)]. Begin with the result that 

A : a; i4 n„. (17) 

onto 

A short proof is given in [H §2.2]. For any t £ [0, T], consider a function u which 
satisfies u{x,t) = for all x S S''"^ = dMd- Define g = AxU. Then 



u{x,t)^ G{x,y)g{y,t) dy, x G Bd, 

JBd 

with G the Green's function for the elliptic boundary value problem 

-Av{x)^g{x), xeMd, 

v{x)=0, xeS'^-^. 
For example, in R^, 

G{x,y)^^log J^-^^y^ x,yeM2, 

2tt \T{x) - y\ 

with T(x) the inverse of x with respect to the unit circle §^. Let gn {x,t) be 
the polynomial referenced in the preceding Lemma [31 and define 

Wnix,t)= G{x,y)gniy,t) dy, a; e B^. (18) 

From (|17p . w„ (•, t) € Xn , < t <T; and {y„ is an approximation of the original 
function u. 



Lemma 4 Assume u{-,t) € C'^'T (id) forO<t<T, with fc > 2, < 7 < 1. 
Then for n > 1, the function Wn (x, t) of I118\) is of the form 



fc=l 



(a;,t) = ^afe(t)?^fe(a;) (19) 

and it satisfies 

\\ui-,t)-w^{-,t)\\c(M.)< ^'f^][f c,,,{g), (20) 

\\y,[u{;t)-w,,i-,t)]\\c{M,) < ^'^C^f ^M(g). (21) 

6fc,-yai (G) 

^^^I+^3^^'=.^(5) (22) 



9 ^ 

— [M(-,t)-W„(-,f)] 



/or < t < T. T/ie constants ai and 012 are given by 
ai(G) = max/ \G{x,y)\dy, 

a2 (G) = max / \VxG {x,y)\ dy, 

x<£Bd JMd 

and these are easily shown to be finite. The remaining constants bk,-f and Ck,-f (<?) 
are taken from Lemma\^ 

Proof. For the error in approximating u, we have 

u (x, t) - Wn {x, t) ^ I G (x, y) [g (y, t) ~ g^ (y, t)] dy, 

JBd 

Va; [u (•, t) - w„ (•, t)] = / Va;G {x, y) [g {y, t) - g^ {y, t)] dy, 

JBd 

— [w (•, t) - Wn (■, 0] = / G {x, y) — [g {y, t) - 5„ {y, t)] dy 
at JBd ot 

Thus 

\\u{-,t) - Wn (-,0110(1^) ^ "1 iG) \\g{-,t)-gn (•,i)llc(i^) 

showing (gni); and ^FQ and ^F^ fohow similarly. ■ 

These results can be extended to the approximation of m (•, i) over J7, by the 
subspace Xn- 

Lemma 5 Assume u{-,t) € G^'^ (H) forO<t<T, with fc > 2, < 7 < 1; 
and assume $ G G"* (B^) with m> fc + 3. Then for n> 1 there exists 






(s,i)=^afc(i)^fe(s), sen, 0<t<T, (23) 



for which 



l|u(-,i)- Wn(-,Ollc(n) 



7T\ < 



\\V^[ui-,t)-Wn{-,t)] 

d 

— [u{-,t)-Wn{,-,t)] 



llc(n) ^ 



< 



c{n) 



UJ2{k,J,u) 

uj3{k,j,u) 

„fe+7-2 



(24) 
(25) 
(26) 



forO<t<T. 



Proof. Use the transformation s = $ (a;) to move between functions over 
fl and functions over B^. By means Lemma [T] for the transformation $, these 
resuhs follow immediately from Lemma |4l ■ 

Combining these results with the Douglas and Dupont Theorem [2] leads to 
the following convergence result for the Galerkin method P^ - ([T5|l . 

Theorem 6 Assume that the solution u of the parabolic problem (l$-(^ satis- 
fies u{-,t) £ C^''' (H) forO<t<T, with k>2, <"f <1. Moreover, assume 
the transformation $ G C™ (B^j) with m > k + 3. Then for n > 1, the solution 
Un of \13 \) - n^) satisfies 



-'n|lL2xL°°' 



'•nwrn^L^ 



= O ( ^-^^-+^-2) 



3 Implementation issues 

Recall the method ([T2|) - (|15p and the notation used there. For notation, let 

aw {t) = [ai {t) ,. .. ,aN (i)] . 
The system (fT51) can be written symbolically as 

Gn^N (t) = Bn {t, Un) 3^ [t) + f^ {t, U„) , 



N 



G„ = [{'4'k,'4'i)]k,i=i^ 



d „ 

{Bn{t,Un))ki^ - y2 / aij{s,t,Unis,t)) 

■ 1 in 

/at {t, ■Un)g = (/ (•, i, U„ {■,t)) , tp() , 



dipk{s,t) dij)i>{s,t) 



dsi 

i-- 



ds 



ds. 



1, 



,N. 



For the implementation, we discuss separately the cases of il C 
ri C M^. In both cases we must address the following issues 

Al. Select a basis {■01, . . . , ^/jjv} for Xn- 

A2. Discuss the numerical integration of the integrals in (P5)) - (PII)) . 



(27) 
(28) 

(29) 

(30) 

|2 and 



A3. Approximate the initial value uq by some uo,n € A'n, as suggested in (IT?)) . 
A4. Discuss the solution of the nonlinear system of differential equations ^7} . 
A5. Evaluate the solution u„ at points of O for each given t. 

Several of these issues were addressed in the previous papers [3] , [5] , [7] , and 
we refer to the discussion in those papers for more complete discussions. 

3.1 Two dimensions 

Let n„ (B2) denote the restriction to 182 of the polynomials over M^. To con- 
struct a basis for the approximation space Xn of pO)) . begin by choosing an 
orthonormal basis {^pi, ■ ■ ■ , ^n} for n„ (B2), using the standard inner product 
for L2 (B2). The dimension of n„ (B2) is 

N = Nn = ^in+l){n + 2) 

There are many possible choices of an orthonormal basis, a number of which are 
enumerated in [T3} §2.3.2] and [321 §l-2]- We have chosen one that is particularly 
convenient for our computations. These are the 'ridge polynomials' introduced 
by Logan and Shepp [12] for solving an image reconstruction problem. We 
summarize here the results needed for our work. 
Let 

v„ = {Pe n„ (B2) : (P, g) = VQ e n„_i} 

the polynomials of degree n that are orthogonal to all elements of n„_i (B2). 
Then the dimension of Vn is n + 1; moreover, 

n„(B2) = Vo®Vi®---®V„ (31) 

It is standard to construct orthonormal bases of each Vn and to then combine 
them to form an orthonormal basis of n„ (B2) using the latter decomposition. 
As an orthonormal basis of V„ we use 

_ 1 n 

<fn.k{x) = —;=Un(xi cos (kh) + X2 sin (kh)) , X e D, h= (32) 

^yn n+1 

for fc = 0, 1, . . . , n. The function U„ is the Chebyshev polynomial of the second 
kind of degree n: 

^^ sm^n+l^^ t^cosff, -1 < t < 1, 71 = 0,1,... 

smO 

The family {'Pn,k}2^Q is an orthonormal basis of V„. As a basis of n„, we order 
{^n.k} lexicographically based on the ordering in (|32|) and (|3T|) : 

{'Pe}e=l = {^0,0, ^1,0, 'Pi,!, 'P2.,0, ■ • ■ , <?n,0, ■ • ■ , ^n,n} 



Returning to ([TU| . we define 

and the basis {V'm.fc : < fc < m, < m < n} for <%"„ is defined using ([TU|. 

V'm.fc (s) = '0n,fc(a;), s==$(a:). 

We will also refer to this basis as {-01, • • • jV'iv}- In general, this is not an 
orthonormal basis; but the hope is that {^i}^^i being orthonornial will result 
in a reasonably well-conditioned matrix for the linear systems associated with 
the solution of P^ . Examples of this for elliptic problems are given in [5], [5], 

To calculate the first order partial derivatives of ipn,k{x), we need U!^{t). 
The values of Un{t) and U^{t) are evaluated using the standard triple recursion 
relations 

Un+l{t)=2tUn{t)-U.a-l{t) 

Ui+,{t) = 2C/„(t) + 2<(t) - Ui_,{t) 

Second derivatives, if needed, can be evaluated similarly. 

For the integrals in ((T5)) . for any dimension d > 2, we first transform them 
to integrals over B,;. For an arbitrary function g defined on fJ, use the transfor- 
mation s = $ (x) to write 

g{s)ds^ / g {^ (x)) det J (x) dx 
with J (x) the Jacobian matrix ([7]) for $(a;). Applying this to the integrals in 



{'ipk,ipi)= / -ijjk (s) ipe (s) ds = / ipk (x) ipe (x) det J (x) dx (33) 
Jn JMd 

( ■^ ( ■^'^^^^k{t)ipk ] ,ipi] 

= / / $ (a;) , i, ^afc (t) V-fc (x) 0fe (a^) det J (x) dx 
•'^■^ \ fc=i / 

> / flij (s,<, u„(s,i))— ^ ^ ds 



iT 



{V0fe(s)}' A(s,w„(s,i)){V^^(s)}ds 



/ |Vi/'fc(a:)} a| a;,i,^afe(t)Vfe(a:;) j|vV'^(a;)}det J(a;) dx 

(35) 



10 



with ^ 

A (a;, t,z)^J (x)"^ A ($ (x) , t, z) J (x)"^ . (36) 

For the numerical approximation of the integrals in psp -ps p with d — 2, 
the integrals being evaluated over the unit disk B2, write a general function g 
as 

g {x) —g{r,9) = g {r cos 6, r sin 0) . 

Then use the formula 

with q > 1 an integer. Here the numbers uii are the weights of the (g + l)-point 
Gauss-Legendre quadrature formula on [0, 1]. The formula ((37| uses the trape- 
zoidal rule with 2q+l subdivisions for the integration over B2 in the azimuthal 
variable. This quadrature ([37l) is exact for all polynomials g € Il2q (IB2). 

To approximate the initial condition uq, as in (jl4[) . we approximate uq ($ (x)) 
by its orthogonal projection onto A'„, 

Vn (mo ° $) = ^ Pj^J 
J = l 

The coefficients {/3j} are obtained by solving the linear system 

Y^ Pj (V'j , V-^) = (wo o $, V^,;) , i = 1, . . . , Af„. (38) 

i=i 

We approximate further by applying the numerical integration l\H\ to each of 
the inner products in this system. With q > n + 2, the matrix coefficients for 
the left side of this linear system will be evaluated exactly. The result of solving 
this system with the associated numerical integration yields an approximation 
to Mo (^ (2;)); and using s — ^ (x), we have an initial estimate of the form given 
in (fT4| . 

To solve the system of ordinary differential equations ([13]), we have used 
the Matlab program odel5s, which is based on the multistep BDF methods 
of orders 1 through 5; see [3, §8.2], [211 P- 60]. In general, there is often 
stiffness when solving differential equations that arise from using a method of 
lines approximation for parabolic problems, and that is our reasoning for using 
the stiff ode code ode 15s rather than an ordinary Runge-Kutta or multistep 
code. No difficulty arose in solving any of our examples when using this code, 
although further work is needed to know whether or not a stiff ode code is 
indeed needed. In our numerical examples, we will give some data on condition 
numbers that arise in our method. 
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3.2 Three dimensions 

Here we denote by n„(B3) the restriction to B3 of polynomials over M.^ of degree 
n or less. The first difference to the two dimensional case is that the dimension 
of n„(B3) is given by 

N = N„ = -{n + l)(n + 2)(n + 3). 

But as with the two dimensional case, there is a wide range of orthonormal basis 
fmictions; see \)A\. We choose the following orthormal basis for n„(B3) 

~ / \ (0-'"^2j+i), I |2 1 N ri / \ 



^771. . 7 U^ 



^_2j „(0,™-2j+i) 



X 



•™,, ,-, pI'"-^'^--'{2\x\' - 1)5^,™_2, (^^ j (39) 

J = 0, . . . , [to/2J , /3 = 0, 1, . . . , 2(m - 2j), m = 0, . . . , n 

The constants Cm,j = 24+t-^J normalize the functions to length one. The 
functions p('''™~2J+2) are the normalized Jacobi polynomials on the interval 
[—1,1] with respect to the inner product 

{V,W) = f (1 +t)"-2j+i y(^t)w{t) dt 

Finally the functions 5'^,m-2j are spherical harmonic functions given by 

{cos (§(J))t^ {cose), P even, 
sin f ^(/.j Tj^^ ^ (cose), /3odd 

Here the constant c^.fc is chosen in such a way that the functions are orthonormal 
on the unit sphere S^ in M^, 



52 



Sp,k{x)S^~^{x)dx = 5p~p5^~^. 



The functions T^, arc the associated Legendre polynomials; see [T8|. In [15], [27] . 
one can also find recurrence formulas for the numerical evaluation of Jacobi and 
Legendre polynomials and their derivatives. 

The bases for the spaces Xn and Xn defined in ^ and (ITUl) are again, see 
dg and dlOl), defined by 

'ipra,],k{x) = (1 - \x\^) ^rn.jM^) (^0) 

■0mj,fc(s) = ^mj\fc(a;), s = >I>(a;) (41) 

For the numerical implementation we can also order the bases in lexicographical 
order (still using the notation "0 and "0), so in the following we can assume that 
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Figure 1: The region fl associated with (|46| and the mapping $ 



we have bases {ipi \ I = 1, ■ • ■ , Nn} and {ipi \ I — 1, . . . , Nn} of Xn and <%"„. 
AU integrals which arise in the formulas (f27|) ~(p0 l) for the approximate solution 
of ([T3l) are transformed to B3 as has been done in (p3l) -(|35 |) . To evaluate the 
resulting integrals over the unit ball in M.^ we use spherical coordinates, and a 
quadrature formula Qq 

g{x)dx= / / / g{r,6,(j))r sni{(f))d(l) dO dr 
,3 Jo Jo Jo 

^ Qqia], where 

2q q q /r* + 1 

QJ5] = 51 XI 51 -^J'^kg I '' , -z-h arccos(^j) 
1=1 j=i fc=i ^ \ ^ ^'^ 

Here g is the representation of g in spherical coordinates. The quadrature 
formula Qq uses a trapezoidal rule in the 6 direction and weighted Gauss- 
Legendre quadrature formulas in the 4> (weights lUj and nodes arccos(^j)) and r 
direction (weights Vk and nodes (^fe + l)/2), as described in [S]. With the help 
of this quadrature formula we can also define the numerical approximation of 
Mo, see dUl) and HSl), by formula (|38|) . 
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4 Numerical examples 

We begin with planar examples, followed by some problems on regions J7 in R'^. 
The examples will all be for the equation 

du is t) 

Q^ ' ^Au{s,t) + f{s,t,u{s,tj), sen, t>0. (42) 

To help in constructing our examples, we use 

f{s,t,z)^his,t,z) + f2{s,t). (43) 

We choose various /i to explore the effects of changes in the type of nonlinearity; 
and /2 is then defined to make the equation (j42p valid for any given u, 

f2{s,t) = ^^^^-{Au{s,t) + h{s,t,u{s,t))},sen, t>0. (44) 



In the reformulation (|35|) . A = I and thus 

A{x,t,z)^ J{x)~'^J{xy'^ . (45) 

4.1 Planar examples 

Begin with the region 57 whose boundary is a limacon. In particular, consider 
the boundary 

(p{0) ^ pi0) (cos 0, sine), , . 

p{9) = 3 + cos9 + 2sm9, 0<9<2tt. ^ ' 

Using the methods of [6], we obtain a mapping $ : B2 -^- il. Each component 
of $ is a polynomial of degree 3. To illustrate the mapping we show the images 
in fl of uniformly spaced circles and radial lines in B2 ; sec Figure [T] and note 
that ri is almost convex. 

As a particular example for solving (|42l) . let 



/i (s, t,z)^e '- cos (vri) , (47) 

u (s, t) = {\-x\- xl) cos {t + O.OSttsiSs) (48) 

with s = $ (x) . For the numerical integration in p7p , q = 2n was chosen, 
where n + 2 is the degree of the approximation Un . This choice of q has always 
been more than adequate, and a smaller choice would often have sufhced. 

To have a time interval of reasonable length, the problem was solved over 
< i < 20, although something longer could have been chosen as well. The 
error was checked at 801 points of ft, chosen as the images under $ of 801 
points distributed over 182- The graph of U12 (-,20) is given in Figure [U and 
the associated error is given in Figure[3l in addition, ||m (-,20) — U12 (-,20)11^ = 
1.94:E — 4. Figure |4] shows the error norm ||u(-,i) — U12 (•,i)||oc for 200 evenly 
spaced values of t in [0, 20]. There is an oscillatory behaviour which is in keeping 



14 




Figure 2: The approximating solution ui2 (s, 20) for the true solution u {s, 20) 
of EHl) over n 



with that of the solution u. To illustrate the spectral rate of convergence of the 
method, Figure [5] gives the error as the degree n varies from 6 to 20. The linear 
behaviour of this semi-log graph implies an exponential rate of convergence of 
Un to u as a function of n. 

An important aspect on which we have not yet commented is the condition- 
ing of the matrices in the system ([27l) . In our use of the Matlab program 
odel5s, we have written (1271) in the form 



(49) 



a^(t) = G„^S„(t,M„)aAr (t) + G„VAr (i,w„), 



The matrix G„ ^-B„ (i, «„) is the Jacobian matrix for this system. Investigating 
experimentally, 

cond {G-^B„) = O {Nl) (50) 



where Nn is the number of equations in (j49|) . As support for this assertion. 
Figure [5] shows the graph of log (-/V^) vs. log (cond (G~^B„)). There is a clear 
linear behaviour and the slope is approximately 1, thus supporting ([50)) . When 
VL is the unit disk, and $ = /, the result ([50]) is still valid experimentally. 

As a second example, one for which 17 is much more nonconvex (although 
still star-like), consider the region J7 with the given boundary function 



ip{e) =p (61) (cos 61, sin 61), 

p (6') = 5 + sin 6 + sin 36* — cos ! 



Q<e <2iT. 



(51) 
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Figure 3: The error in the approximating solution ui2 (s, 20) for the true solution 
u (s, 20) of dm) over n 



As before, an extension $ to B2 is constructed using the methods of 6 . The 
mapping $ is a polynomial of degree 7 in each component; and the images in O 
of uniformly spaced circles and radial lines in B2 are shown in Figure [T] 

Again, use the function /i of (|T7)) and the solution u of (|35]). The solution 
U20 (',20) is shown in Figure|8]over this new region, and ||u (-,20) — U20 (•, 20)11^^ 
= 0.00136. Figure [9] shows the error in U20 (', t) over time, and Figure [TOl shows 
how the error in m„ varies with the degree n. The latter again indicates a 
spectral order of convergence, although slower than that shown in Figure El 
The condition numbers still satisfy the empirical estimate of (|50p. 



4.2 A three-dimensional example 

Here we will study one domain Q which we investigated already in a previous 
article for the purpose of analyzing the spectral method for Dirichlet problems; 
see [2]. The domain has the advantage that the transformation $ is known 
throughout B3 and even the inverse transformation ^ is known explicitly. The 
knowledge of ^ is not necessary for the use of the spectral method but makes 
the construction of an explicit solution easier. The mapping $ : B3 i-->- fi, 
(si,S2,S3) = $(xi,X2,a;3) is given by 



Sl = 


= Xi - 2:2 - 


S2 = 


= Xi + X2 


S3 = 


= 2x3 + bx 



(52) 
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Figure 4: The error ||m (•, i) — ui2 (•, t)\\^ for the true solution u (s, t) of ([48] 



where < a, 6 < 1 are two parameters. Figures [TTJ [T^ show an example of the 
surface of fl from two different angles. The inverse ^ : fJ i— > B3 is given by 



i 

Xi ^ - 
a 

1 

X2 = - 

a 
1 


-1 + v/l + a(si + S2) 




as2 + l- v/l + 


a{si + S2) 


-1 + ^/l + bs3 





Furthermore the Jacobian for $ is given by 

/ l + 2axi -1 

J{x) =1 1 

\ 2 + 26x3 

with determinant 

det(J(a;)) = 4(1 + axi){l + bxs). 

This allows us also to calculate A, see (|45p . directly 



/ 



1 



A{x) 



axi 



2(l + axi)2 2(l + aa;i)2 

aa;i 1 + axi + 2a?'x\ 

2{l + axif 2(l + aa;i)2 





4(1 + 5x3)2 / 
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Figure 5: n vs. max llu (•, t) — m„ (•, t)|l 

0<t<20 ^ 



Again we use the spectral method to solve (UU where / is given by P5)) and 
(|44|) . As a particular example for solving (|42p . let 



fi{s,t,z) = e ^cos(7ri), 

u(s, t) = {1- x\~ x\- x\) cos(t + 0.057rsiS2S3) (53) 

where (xi, a;2, 2:3) = ^(si, S2, S3) with a = 0.7 and b = 0.9. Numerical results are 
given in Figures [131 [TH Figure [15] seems to indicate that the relation (|50| for the 
condition number of the Jacobian G~^Bn is also valid in the three dimensional 
case. 
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Figure 9: The error \\u{-,t) — U2o{-,t)\\^ for the true sohition u{s,t) of 
over the il of Figure [7] 
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Figure 10: n VS. max llu (•,i) — u„ (•, i)|l for the il of FigurelT] 

0<t<20 °° 
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Figure 11: The surface of il. with parameters a = 0.7 and b = 0.9; see ([52 
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Figure 12: The surface of J7, see ([52l) . seen from the z-axis 
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Figure 13: The error \\u{-,t) — ui2 (•,i)||oo for the true solution u{s,t) of (|55|) 
over the 57 of Figure [TT] 
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Figure 14: n vs. max |lu(-,i) — m„ (-.illl^ over the fl of Figure [Til 

*^ 0<t<20" ^ ^ ^ ^"°° ^ 
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Figure 15: log (iV^) vs. log (cond {G-^Bn)) for the n of FigureE] 
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